1 C T-00 PROGRAM
2 DIMENSION NDLTH(4),CDH(4),AR(5,5),AI(5,5),BR(5,5),BI(5,5),PN(6),
3 1X(10),Y(10),T(10,10),RE(10,10),IM(10,10),TN(6,10,10)
4 DOUBLE PRECISION THETA,CDH,DLTH,SRMUL,STH,CTH,KR,DKR,PN,FCT,X,Y,BT
5 1,CT,W1,W2,W3,W4,W5,AR,AI,BR,BI,RE,IM,K1,R,DR,ALPHA,A,B,C,D,PI,
6 1THETA1,THETA2,THETA3,THETA4,E,F,G
7 COMPLEX*16 T,TN
8 COMMON/FAC/FCT(100),NRISK,NTRIX
9 NAMELIST/HEJ/AR,AI,BR,BI,RE,IM,T
10
11 open (6, file='matrix1.dat', status='unknown')
12 cc open (11, file='matrix2.dat', form='unformatted', status='unknown')
13
14 cc CALL UNSTAE
15 PI = 4.0D0*DATAN(1.0D0)
16 NRISK = 57
17 NTRIX = 39
18 NEND = 78
19 NR = 5
20 NP = 0
21 ISYM = 0
22 K1 = 0.5D0
23 c diameter
24 A = 1.00D0
25
26 C = -0.1D0
27 NSECT = 1
28 NDLTH(1) = 192
29 THETA1 = PI
30 CDH(1) = THETA1/DFLOAT(NDLTH(1))
31 N = NRISK
32 L = NTRIX
33 M = NEND-2
34 N1 = N-1
35 DO 5 K = 1,100
36 5 FCT(K) = 0.0D0
37 FCT(1) = 1.0D0
38 DO 6 K = 1,N1
39 6 FCT(K+1) = FCT(K)*DFLOAT(K)
40 FCT(N+1) = 0.1D0**L*FCT(N)*DFLOAT(N)
41 DO 7 K = N,M
42 7 FCT(K+2) = FCT(K+1)*DFLOAT(K+1)
43 ND = 2*NR
44 NR1 = NR+1
45 DO 150 M1 = 1,NR1
46 M = M1-1
47 MD = M
48 IF(M.EQ.0) MD = 1
49 DO 25 J = 1,NR
50 DO 25 I = 1,NR
51 AR(I,J) = 0.0D0
52 AI(I,J) = 0.0D0
53 BR(I,J) = 0.0D0
54 BI(I,J) = 0.0D0
55 25 CONTINUE
56 THETA = 0.0D0
57 DO 90 ISECT = 1,NSECT
58 NTHETA = NDLTH(ISECT)+1
59 DLTH = CDH(ISECT)
60 NFIX = 1
61 DO 89 K = 1,NTHETA
62 IF(K-1)31,31,32
63 31 CONTINUE
64 SRMUL = DLTH*7.0D0/22.5D0
65 IF(ISECT.EQ.1) GO TO 89
66 C Seite 1
67
68
69 GO TO 41
70 32 CONTINUE
71 IF(K-NTHETA)34,33,33
72 33 CONTINUE
73 SRMUL = DLTH*7.0D0/22.5D0
74 IF(ISECT-NSECT)40,89,89
75 34 CONTINUE
76 GO TO (35,36,37,38),NFIX
77 35 CONTINUE
78 SRMUL = DLTH*32.0D0/22.5D0
79 NFIX = 2
80
81
82 GO TO 39
83 36 CONTINUE
84 SRMUL = DLTH*12.0D0/22.5D0
85 NFIX = 3
86 GO TO 39
87 37 CONTINUE
88 SRMUL = DLTH*32.0D0/22.5D0
89 NFIX = 4
90 GO TO 39
91 38 CONTINUE
92 SRMUL = DLTH*14.0D0/22.5D0
93 NFIX = 1
94 39 CONTINUE
95 40 CONTINUE
96 THETA = THETA+DLTH
97 STH = DSIN(THETA)
98 CTH = DCOS(THETA)
99 CALL LEG(THETA,M,NR1,PN)
100 41 CONTINUE
101 CALL TRCIR(STH,CTH,C,A,R,DR)
102 KR = K1*R
103 DKR = K1*DR
104 IF(K.EQ.1) GO TO 52
105 CALL BN(KR,NR1,X,Y)
106 BT= DABS(KR*KR*(X(2)*Y(1)-X(1)*Y(2))-1.0D0)
107 CT= DABS(KR*KR*(X(NR1)*Y(NR)-X(NR)*Y(NR1))-1.0D0)
108 IF(BT-1.D-10)48,48,45
109 45 CONTINUE
110 PRINT 46
111 46 FORMAT ('0 ERROR IN BESSEL TEST')
112 PRINT 47,BT,ISECT,K
113 47 FORMAT(D25.15,2I5)
114 48 CONTINUE
115 IF(CT-1.D-10)52,52,49
116 49 CONTINUE
117 PRINT 50
118 50 FORMAT ('0 ERROR IN NEUMAN TEST')
119 PRINT 51,CT,ISECT,K
120 51 FORMAT(D25.15,2I5)
121 52 CONTINUE
122 DO 88 I = MD,NR
123 DO 88 J = MD,NR
124 I1 = I+1
125 J1 = J+1
126 IF(M.EQ.0) GO TO 74
127 IF(ISYM.GT.0.AND.MOD(I+J,2).EQ.0) GO TO 74
128 C Seite 2
129
130
131 W1 = DFLOAT(I+J)*CTH*PN(I1)*PN(J1)-DFLOAT(I+M)*PN(I)*PN(J1)-DFLOAT
132 1(J+M)*PN(I1)*PN(J)
133 W2 = KR*KR*X(I1)
134
135
136 W3 = W1*W2
137 IF(ISYM.NE.2) GO TO 71
138 IF(J-1)72,71,71
139 71 CONTINUE
140 BI(I,J) = BI(I,J)+SRMUL*W3*Y(J1)/STH
141 72 CONTINUE
142 IF(J-I)74,73,73
143 73 CONTINUE
144 BR(I,J) = BR(I,J)+SRMUL*W3*X(J1)/STH
145 74 CONTINUE
146 IF(ISYM.GT.0.AND.MOD(I+J,2).NE.0) GO TO 84
147 W4 = KR*(KR*X(I)-DFLOAT(I)*X(I1))
148 W5 = DFLOAT(I1)*DKR*STH*X(I1)
149 W1 = PN(I1)*PN(J1)*(W4*(DFLOAT(I*J)*(CTH**2)+DFLOAT(M*M))+DFLOAT
150 1(I*J)*W5*CTH)
151 W2 = DFLOAT(I*(J+M))*PN(I1)*PN(J)*(CTH*W4+W5)
152 W3 = DFLOAT(I*M)*PN(I)*W4*(DFLOAT(J+M)*PN(J)-DFLOAT(J)*CTH*PN(J1))
153 IF(ISYM.NE.2) GO TO 81
154 IF(J-I)82,81,81
155 81 CONTINUE
156 AI(I,J) = AI(I,J)+SRMUL*(W1-W2+W3)*Y(J1)/STH
157
158
159
160 82 CONTINUE
161 IF(J-I)84,83,83
162 83 CONTINUE
163 AR(I,J) = AR(I,J)+SRMUL*(W1-W2+W3)*X(J1)/STH
164 84 CONTINUE
165 88 CONTINUE
166 89 CONTINUE
167 90 CONTINUE
168 DO 92 I = 2,NR
169 JEND = I-1
170 DO 91 J = 1,JEND
171 BR(I,J) = BR(J,I)
172 AR(I,J) = AR(J,I)
173 IF(ISYM.NE.2) GO TO 91
174 BI(I,J) = BI(J,I)
175 AI(I,J) = AI(J,I)
176 91 CONTINUE
177 92 CONTINUE
178 DO 97 I = MD,NR
179 DO 97 J = MD,NR
180 I1 = I+1
181 J1 = J+1
182 W1 = -DSQRT(DFLOAT((2*I+1)*(2*J+1))/DFLOAT(I*I1*J*J1))/2.0D0
183 IF(M)96,96,95
184 95 CONTINUE
185 W1 = W1*DSQRT(FCT(I1-M)*FCT(J1-M)/(FCT(I1+M)*FCT(J1+M)))
186 96 CONTINUE
187 AR(I,J) = W1*AR(I,J)
188 AI(I,J) = W1*AI(I,J)
189 BR(I,J) = DFLOAT(M)*W1*BR(I,J)
190 BI(I,J) = DFLOAT(M)*W1*BI(I,J)
191 97 CONTINUE
192 CALL PERFT(NP,M,NR,ND,AR,AI,BR,BI,T,RE,IM,X,Y)
193 C Seite 3
194
195
196 DO 130 I = 1,ND
197 DO 130 J = 1,ND
198 TN(M1,I,J) = T(I,J)
199 130 CONTINUE
200 IF(M.NE.1) GO TO 150
201 WRITE(6,HEJ)
202
203
204 150 CONTINUE
205 WRITE(11) TN
206 PRINT 161
207 161 FORMAT ('0 TN NOW HOPEFULLY REPLACED INTO DATA SET')
208 STOP
209 cc DEBUG SUBCHK
210 END
211 C Seite 4
212
213
214
215
216 SUBROUTINE BESSEL(NORDER,ARGMNT,ANSWR,IERROR)
217 DOUBLE PRECISION ARGMNT,ANSWR,X,SUM,APR,TOPR,CI,CNI,ACR,PROD,FACT
218 IERROR = 0
219 N = NORDER
220 X = ARGMNT
221 SUM = 1.0D0
222 APR = 1.0D0
223 TOPR = -0.5D0*X*X
224 CI = 1.0D0
225 CNI = DFLOAT(2*N+3)
226 DO 60 I = 1,100
227 ACR = TOPR*APR/(CI*CNI)
228 SUM = SUM+ACR
229 IF(DABS(ACR/SUM)-1.0D-20)100,100,40
230 40 APR = ACR
231 CI = CI+1.0D0
232 CNI = CNI+2.0D0
233 60 CONTINUE
234 IERROR = I
235 PRINT 10
236 10 FORMAT(24H0 ERROR IN SUM OF BESSEL)
237 GO TO 200
238 100 PROD = DFLOAT(2*N+1)
239 FACT = 1.0D0
240 IF(N)160,160,120
241 120 DO 140 IFCT = 1,N
242 FACT = FACT*X/PROD
243 PROD = PROD-2.0D0
244 140 CONTINUE
245 160 ANSWR = FACT*SUM
246 200 RETURN
247 END
248 C Seite 23
249
250
251 SUBROUTINE BN(PCKR,NRINK,BSSLSP,CNEUMN)
252 DIMENSION BSSLSP(NRINK),CNEUMN(NRINK)
253 DOUBLE PRECISION PCKR,ANSWR,ANSA,ANSB,ANSC,CONN,CMULN,SNSA,SNSB,
254 1SNSC,BSSLSP,CNEUMN
255 NRANKI = NRINK
256 NRANK = NRANKI-1
257 NVAL = NRANK-1
258 DO 40 I = 1,4
259 CALL BESSEL(NVAL,PCKR,ANSWR,IERROR)
260 IF(IERROR)20,20,32
261 20 ANSA = ANSWR
262 NVAL = NVAL+1
263 CALL BESSEL(NVAL,PCKR,ANSWR,IERROR)
264 IF(IERROR)24,24,28
265 24 ANSB = ANSWR
266 GO TO 60
267 28 NVAL = NVAL-1
268 32 NVAL = NVAL+NRANK
269 40 CONTINUE
270 STOP
271 60 IF(NVAL-NRANK)100,100,64
272 64 IEND = NVAL-NRANK
273 CONN = DFLOAT(2*(NVAL-1)+1)
274 DO 72 IP = 1,IEND
275 ANSC = CONN*ANSA/PCKR-ANSB
276 CONN = CONN-2.0D0
277 ANSB = ANSA
278 ANSA = ANSC
279 72 CONTINUE
280 100 BSSLSP(NRANKI) = ANSB
281 BSSLSP(NRANKI-1) = ANSA
282 CONN = DFLOAT(NRANK+NRANK-1)
283 IE = NRANKI-2
284 JE = IE
285 DO 120 JB = 1,JE
286 ANSC = CONN*ANSA/PCKR-ANSB
287 BSSLSP(IE) = ANSC
288 ANSB = ANSA
289 ANSA = ANSC
290 IE = IE-1
291 CONN = CONN-2.0D0
292 120 CONTINUE
293 CMULN = 3.0D0
294 SNSA = -DCOS(PCKR)/PCKR
295 SNSB = -DCOS(PCKR)/(PCKR*PCKR)-DSIN(PCKR)/PCKR
296 CNEUMN(1) = SNSA
297 CNEUMN(2) = SNSB
298 DO 280 I = 3,NRANKI
299 SNSC = CMULN*SNSB/PCKR-SNSA
300 CNEUMN(I) = SNSC
301 SNSA = SNSB
302 SNSB = SNSC
303 CMULN = CMULN+2.0D0
304 280 CONTINUE
305 RETURN
306 END
307 C Seite 24
308
309
310 SUBROUTINE CBESS (NORDER,ARGMNT,ANSWR,IERROR)
311 COMPLEX*16 ARGMNT,ANSWR,X,SUM,APR,TOPR,CI,CNI,ACR,PROD,FACT
312 IERROR = 0
313 N = NORDER
314 X = ARGMNT
315 SUM = (1.0D0,0.0D0)
316 APR = (1.0D0,0.0D0)
317 TOPR = -(0.5D0,0.0D0)*X*X
318 CI = (1.0D0,0.0D0)
319 CNI = DCMPLX(DFLOAT(2*N+3),0.0D0)
320 DO 60 I = 1,100
321 ACR = TOPR*APR/(CI*CNI)
322 SUM = SUM+ACR
323 IF(CDABS(ACR/SUM)-1.0D-20)100,100,40
324 40 APR = ACR
325 CI = CI+(1.0D0,0.0D0)
326 CNI = CNI+(2.0D0,0.0D0)
327 60 CONTINUE
328 IERROR = 1
329 PRINT 10
330 10 FORMAT('0 ERROR IN SUM OF CBESS')
331 GO TO 200
332 100 PROD = DCMPLX(DFLOAT(2*N+1),0.0D0)
333 FACT = (1.0D0,0.0D0)
334 IF(N)160,160,120
335 120 DO 140 IFCT = 1,N
336 FACT = FACT*X/PROD
337 PROD = PROD-(2.0D0,0.0D0)
338 140 CONTINUE
339 160 ANSWR = FACT*SUM
340 200 RETURN
341 END
342 C Seite 25
343
344
345 SUBROUTINE CBN(PCKR,NRINK,BSSLSP,CNEUMN)
346 DIMENSION BSSLSP(NRINK),CNEUMN(NRINK)
347 COMPLEX*16 PCKR,ANSWR,ANSA,ANSB,ANSC,CONN,CMULN,SNSA,SNSB,SNSC,
348 1BSSLSP,CNEUMN
349 NRANKI = NRINK
350 NRANK = NRANKI-1
351 NVAL = NRANK-1
352 DO 40 I = 1,4
353 CALL CBESS(NVAL,PCKR,ANSWR,IERROR)
354 IF(IERROR)20,20,32
355 20 ANSA = ANSWR
356 NVAL = NVAL+1
357 CALL CBESS(NVAL,PCKR,ANSWR,IERROR)
358 IF(IERROR)24,24,28
359 24 ANSB = ANSWR
360 GO TO 60
361 28 NVAL = NVAL-1
362 32 NVAL = NVAL+NRANK
363 40 CONTINUE
364 STOP
365 60 IF(NVAL-NRANK)100,100,64
366 64 IEND = NVAL-NRANK
367 CONN = DCMPLX(DFLOAT(2*(NVAL-1)+1),0.0D0)
368 DO 72 IP = 1,IEND
369 ANSC = CONN*ANSA/PCKR-ANSB
370 CONN = CONN-(2.0D0,0.0D0)
371 ANSB = ANSA
372 ANSA = ANSC
373 72 CONTINUE
374 100 BSSLSP(NRANKI) = ANSB
375 BSSLSP(NRANKI-1) = ANSA
376 CONN = DCMPLX(DFLOAT(NRANK+NRANK-1),0.0D0)
377 IE = NRANKI-2
378 JE = IE
379 DO 120 JB = 1,JE
380 ANSC = CONN*ANSA/PCKR-ANSB
381 BSSLSP(IE) = ANSC
382 ANSB = ANSA
383 ANSA = ANSC
384 IE = IE-1
385 CONN = CONN-(2.0D0,0.0D0)
386 120 CONTINUE
387 CMULN = (3.0D0,0.0D0)
388 SNSA = -CDCOS(PCKR)/PCKR
389 SNSB = -CDCOS(PCKR)/(PCKR*PCKR)-CDSIN(PCKR)/PCKR
390 CNEUMN(1) = SNSA
391 CNEUMN(2) = SNSB
392 DO 280 I = 3,NRANKI
393 SNSC = CMULN*SNSB/PCKR-SNSA
394 CNEUMN(I) = SNSC
395 SNSA = SNSB
396 SNSB = SNSC
397 CMULN = CMULN+(2.0D0,0.0D0)
398 280 CONTINUE
399 RETURN
400 END
401 C Seite 26
402
403
404 SUBROUTINE LEG(THETA,M,NRJNK,PNMLLG)
405 DIMENSION PNMLLG(NRJNK)
406 DOUBLE PRECISION THETA,PLA,PLB,PLC,DTWM,CNM,CNMM,CNMUL,PRODM,
407 1QUANM,PNMLLG
408 NRANKI = NRJNK
409 KMV = M
410 DTWM = DFLOAT (2*M+1)
411 IF(THETA)16,4,16
412 4 IF(KMV)12,12,6
413 6 DO 8 ILG = 1,NRANKI
414 PNMLLG(ILG) = 0.0D0
415 8 CONTINUE
416 GO TO 88
417 12 PNMLLG(1) = 1.0D0
418 PLA = 1.0D0
419 GO TO 48
420 16 IF(KMV)20,20,40
421 20 PLA = 1.0D0
422 PLB = DCOS(THETA)*PLA
423 PNMLLG(1) = PLA
424 PNMLLG(2) = PLB
425 IBEG = 3
426 GO TO 60
427 40 DO 44 ILG = 1,KMV
428 PNMLLG(ILG) = 0.0D0
429 44 CONTINUE
430 PRODM = 1.0D0
431 QUANM = DFLOAT(KMV)
432 DO 52 IFCT = 1,KMV
433 QUANM = QUANM+1.0D0
434 PRODM = QUANM*PRODM/2.0D0
435 52 CONTINUE
436 PLA = PRODM*DSIN(THETA)**KMV
437 PNMLLG(KMV+1) = PLA
438 48 PLB = DTWM*DCOS(THETA)*PLA
439 PNMLLG(KMV+2) = PLB
440 IBEG = KMV+3
441 60 CNMUL = DFLOAT(IBEG+IBEG-3)
442 CNM = 2.0D0
443 CNMM = DTWM
444 DO 80 ILGR = IBEG,NRANKI
445 PLC = (CNMUL*DCOS(THETA)*PLB-CNMM*PLA)/CNM
446 PNMLLG(ILGR) = PLC
447 PLA = PLB
448 PLB = PLC
449 CNMUL = CNMUL+2.0D0
450 CNM = CNM+1.0D0
451 CNMM = CNMM+1.0D0
452 80 CONTINUE
453 88 RETURN
454 END
455 C Seite 27
456
457
458 FUNCTION TRIXJ(J1,J2,J,M,FCT,N,L)
459 IMPLICIT REAL*8(A-H,O-Z)
460 DIMENSION FCT(1)
461 INTEGER Z,ZMIN,ZMAX
462 Y=1.0D0
463 CC=0.0D0
464 JSUM=J1+J2+J
465 JM1=J1-IABS(M)
466 JM2=J2-IABS(M)
467 IF((MOD(JSUM,2).NE.0).OR.(MOD(JM1,2).NE.0).OR.(MOD(JM2,2).NE.0)
468 1.OR.(JM1.LT.0).OR.(JM2.LT.0))
469 2GO TO 3
470 IF((J.GT.J1+J2).OR.(L.LT.IABS(J1-J2))) GO TO 1
471 ZMIN=0
472 IF(J-J2+M.LT.0) ZMIN=-J+J2-M
473 IF(J-J1+M+ZMIN.LT.0) ZMIN=-J+J1-M
474 ZMAX=J1+J2-J
475 IF(J2-M-ZMAX.LT.0) ZMAX=J2-M
476 IF(J1-M-ZMAX.LT.0) ZMAX=J1-M
477 JA=(J1+M)/2+1
478 JB=JA-M
479 JC=(J2-M)/2+1
480 JD=JC+M
481 JE=J/2+1
482 JF=(J1+J2-J)/2+1
483 JG=JA+JB-JF
484 JH=JC+JD-JF
485 JJ=2*JE+JF-1
486 IF(JJ.GT.N) Y=0.1D0**L
487 IF(FCT(JJ))7,5,7
488 7 CONTINUE
489 IA=ZMIN/2
490 IB=JF-IA+1
491 IC=JB-IA+1
492 ID=JC-IA+1
493 IE=JA-JF+IA
494 IF=JD-JF+IA
495 FASE=1.0D0
496 IF(MOD(IA,2).EQ.0) FASE=-FASE
497 Z=ZMIN
498 2 IA=IA+1
499 IB=IB-1
500 IC=IC-1
501 ID=ID-1
502 IE=IE+1
503 IF=IF+1
504 FASE=-FASE
505 CC=CC+FASE/(FCT(IA)*FCT(IB)*FCT(IC)*FCT(ID)*FCT(IE)*FCT(IF))
506 Z=Z+2
507 IF(Z.LE.ZMAX) GO TO 2
508 FASE=-DSIGN(1.0D0,CC)
509 IF(MOD(J1-J2,4).EQ.0) FASE=-FASE
510 CC=FASE*DSQRT(Y*FCT(JB)*FCT(JC)*FCT(JE)*CC*FCT(JA)*FCT(JD)*FCT(JE)
511 1*CC*FCT(JF)*FCT(JG)*FCT(JH)/FCT(JJ))
512 1 TRIXJ=CC
513 RETURN
514 3 PRINT 4
515 4 FORMAT('0 ERROR IN ARGUMENT OF TRIXJ')
516 CALL EXIT
517 5 PRINT 6
518
519
520 6 FORMAT('0 ERROR FACTORIALS EXCEEDED IN TIXJ')
521 CALL EXIT
522 END
523 C Seite 28
524
525
526 SUBROUTINE VPSI(RAD,THETA,FHI,NP,NB,M,ND,PSIR,PSITH,PSIFH,PN,U,V)
527 DIMENSION PSIR(1),PSITH(1),PSIFH(1),PN(1),U(1),V(1)
528 DOUBLE PRECISION RAD,R,THETA,T,FHI,F,PN,U,V,FCT,FR,FR1,FR2,B,C
529 COMPLEX*16 FC,FC1,FC2,PSIR,PSITH,PSIFH
530 COMMON/FAC/FCT(100),NRISK,NTRIX
531 R = RAD
532 T = THETA
533 F = FHI
534 NP1 = NP-1
535 K = M
536 N = ND
537 L = N/2
538 L1 = L+1
539 L2 = L+2
540 FC = (0.0D0,1.0D0)
541 IF(NB.EQ.1) GO TO 5
542 CALL BN(R,L1,U,V)
543 B = DABS(R*R*(U(2)*V(1)-U(1)*V(2))-1.0D0)
544 C = DABS(R*R*(U(L1)*V(L)-U(L)*V(L1))-1.0D0)
545 PRINT 3
546 3 FORMAT('0 BESSEL- NEUMAN- TEST FOR VPSI')
547 PRINT 4,B,C
548 4 FORMAT(2D25.15)
549 5 CONTINUE
550 CALL LEG(T,K,L2,PN)
551 DO 42 I = 1,L
552 I1 = I+1
553 I2 = I+2
554 IF(I-K)6,7,7
555 6 FR = 0.0D0
556 GO TO 10
557 7 IF(K)8,8,9
558 8 FR = DSQRT(DFLOAT(2*I+1)/(DFLOAT(I*I1)*16.0D0*DATAN(1.0D0)))
559 GO TO 10
560 9 FR = DSQRT(DFLOAT(2*I+1)*FCT(I1-K)/(DFLOAT(I*I1)*FCT(I1+K)*8.0D0*
561 1DATAN(1.0D0)))
562 10 CONTINUE
563 IF(T)16,16,19
564 16 CONTINUE
565 IF(K-1)18,17,18
566 17 CONTINUE
567 FR1 = DSQRT(DFLOAT(2*I+1)/(32.0D0*DATAN(1.0D0)))
568 FR2 = DSQRT(DFLOAT(2*I+1)/(32.0D0*DATAN(1.0D0)))
569 GO TO 20
570 18 CONTINUE
571 FR1 = 0.0D0
572 FR2 = 0.0D0
573 GO TO 20
574 19 CONTINUE
575 FR1 = FR*PN(I1)*DFLOAT(K)/DSIN(T)
576 FR2 = -FR*(DFLOAT(I1)*DCOS(T)*PN(I1)-DFLOAT(I1-K)*PN(I2))/DSIN(T)
577 20 CONTINUE
578 GO TO (30,31,32),NB
579 30 CONTINUE
580 FC1 = (-FC)**I1
581 FC2 = (-FC)**I
582 GO TO 33
583 31 CONTINUE
584 C Seite 29
585
586
587 FC1 = DCMPLX(U(I1),0.0D0)
588 FC2 = DCMPLX(U(I)-DFLOAT(I)*U(I1)/R,0.0D0)
589
590
591 GO TO 33
592 32 CONTINUE
593 FC1 = DCMPLX(U(I1),V(I1))
594 FC2 = DCMPLX(U(I)-DFLOAT(I)*U(I1)/R,V(I)-DFLOAT(I)*V(I1)/R)
595 33 CONTINUE
596 PSIR(2*I-1) = (0.0D0,0.0D0)
597 IF(NB-1)34,34,35
598 34 CONTINUE
599 PSIR(2*I) = (0.0D0,0.0D0)
600 35 CONTINUE
601 GO TO (36,39),NP1
602 36 CONTINUE
603 IF(NB-1)38,38,37
604 37 CONTINUE
605 PSIR(2*I) = FC1*DCMPLX(DFLOAT(I*I1)*FR*PN(I1)*DSIN(DFLOAT(K)*F)/R,
606 10.0D0)
607 38 CONTINUE
608 PSITH(2*I-1) = -FC1*DCMPLX(FR1*DSIN(DFLOAT(K)*F),0.0D0)
609 PSITH(2*I) = FC2*DCMPLX(FR2*DSIN(DFLOAT(K)*F),0.0D0)
610 PSIFH(2*I-1) = -FC1*DCMPLX(FR2*DCOS(DFLOAT(K)*F),0.0D0)
611 PSIFH(2*I) = FC2*DCMPLX(FR1*DCOS(DFLOAT(K)*F),0.0D0)
612 GO TO 42
613 39 CONTINUE
614 IF(NB-1)41,41,40
615 40 CONTINUE
616 PSIR(2*I) = FC1*DCMPLX(DFLOAT(I*I1)*FR*PN(I1)*DCOS(DFLOAT(K)*F)/R,
617 10.0D0)
618 41 CONTINUE
619 PSITH(2*I-1) = FC1*DCMPLX(FR1*DCOS(DFLOAT(K)*F),0.0D0)
620 PSITH(2*I) = FC2*DCMPLX(FR2*DCOS(DFLOAT(K)*F),0.0D0)
621 PSIFH(2*I-1) = -FC1*DCMPLX(FR2*DSIN(DFLOAT(K)*F),0.0D0)
622 PSIFH(2*I) = -FC2*DCMPLX(FR1*DSIN(DFLOAT(K)*F),0.0D0)
623 42 CONTINUE
624 RETURN
625 END
626 C Seite 30
627
628
629 SUBROUTINE VKOEF (BHETA,NP,M,ND,AP,PN)
630 DIMENSION PN(1),AP(1)
631 DOUBLE PRECISION BHETA,U,DI,PN,FR,FCT
632 COMPLEX*16 FC1,FC2,fc3,AP
633 COMMON/FAC/FCT(100),NRISK,NTRIX
634 N = NP
635 N1 = N+1
636 K = M
637 L = ND/2
638 L2 = L+2
639 U = BHETA
640 DI = DATAN(1.0D0)
641 FC1 = (0.0D0,1.0D0)
642 CALL LEG(U,K,L2,PN)
643 IF(U)10,10,20
644 10 DO 100 I = 1,L
645 IF(K-1)11,12,11
646 11 AP(2*I-1) = (0.0D0,0.0D0)
647 AP(2*I) = (0.0D0,0.0D0)
648 GO TO 100
649 12 FR = DSQRT(8.0D0*DI*DFLOAT(2*I+1))
650 FC2 = DCMPLX(FR,0.0D0)
651 GO TO (13,14),N1
652 13 AP(2*I-1) = -(FC1**I)*FC2
653 AP(2*I) = -(FC1**(I+1))*FC2
654 GO TO 100
655 14 AP(2*I-1) = (FC1**I)*FC2
656 AP(2*I) = (FC1**(I+3))*FC2
657 100 CONTINUE
658 GO TO 300
659 20 DO 200 I = 1,L
660 I1 = I+1
661 I2 = I+2
662 IF(K)21,21,30
663 21 FR = 4.0D0*DSQRT((DI*DFLOAT((2*I+1)*(I+1)))/DFLOAT(I))*(DCOS(U)*
664 1PN(I1)-PN(I2))/DSIN(U)
665 FC2 = DCMPLX(FR,0.0D0)
666 GO TO (22,23),N1
667 22 AP(2*I-1) = (FC1**I)*FC2
668 AP(2*I) = (0.0D0,0.0D0)
669 GO TO 200
670 23 AP(2*I-1) = (0.0D0,0.0D0)
671 AP(2*I) = (FC1**(I+1))*FC2
672 GO TO 200
673 30 IF(I-K)34,31,31
674 31 FR = 4.0D0*DSQRT((2.0D0*DI*DFLOAT(2*I+1)*FCT(I-K+1))/(DFLOAT(I*(I
675 1+1))*FCT(I+K+1)))
676 FC2 = DCMPLX((FR*(DFLOAT(I+1)*DCOS(U)*PN(I1)-DFLOAT(I-K+1)*PN(I2)
677 1))/DSIN(U),0.0D0)
678 FC3 = DCMPLX((DFLOAT(K)*FR*PN(I1))/DSIN(U),0.0D0)
679 GO TO (32,33),N1
680 32 AP(2*I-1) = (FC1**I)*FC2
681 AP(2*I) = -(FC1**(I+1))*FC3
682 GO TO 200
683 33 AP(2*I-1) = (FC1**I)*FC3
684 AP(2*I) = (FC1**(I+1))*FC2
685 GO TO 200
686 34 AP(2*I-1) = (0.0D0,0.0D0)
687 AP(2*I) = (0.0D0,0.0D0)
688 200 CONTINUE
689
690
691 300 RETURN
692 END
693 C Seite 31
694
695
696 SUBROUTINE VR(AR,NP,M,ND,R1,R2,R3,R4,X1,X2,Y1,Y2)
697 DIMENSION R1(ND,1),R2(ND,1),R3(ND,1),R4(ND,1),X1(1),X2(1),Y1(1),Y2
698 1(1)
699 DOUBLE PRECISION AR,U,V,F3J1,F3J2,F3J3,FACT1,FACT2,FACT3,FACT4,
700 1FACT5,FACT6,FACT7,FACT8,FACT9,FACT10,ACR1,ACR2,ACR3,ACR4,X1,X2,
701 1Y1,Y2,B1,B2,CN1,CN2,FCT,trixj
702 COMPLEX*16 R1,R2,R3,R4,W1,W2,W3,W4,W5,W6,W7,W8
703 COMMON/FAC/FCT(100),NRISK,NTRIX
704 NP1 = NP+1
705 N = M
706 NK = ND+1
707 U = AR
708 V = 0.5D0*AR
709 CALL BN(U,NK,X1,Y1)
710 CALL BN(V,NK,X2,Y2)
711 L1 = NK
712 L = NK-1
713 B1 = DABS(U**2*(X1(2)*Y1(1)-X1(1)*Y1(2))-1.0D0)
714 B2 = DABS(V**2*(X1(2)*Y2(1)-X2(1)*Y2(2))-1.0D0)
715 CN1 = DABS(U**2*(X1(L1)*Y1(L)-X1(L)*Y1(L1))-1.0D0)
716 CN2 = DABS(V**2*(X2(L1)*Y2(L)-X2(L)*Y2(L1))-1.0D0)
717 PRINT 89
718 89 FORMAT('0 BESSEL- NEUMAN- TEST FOR VR')
719 PRINT 90, B1,CN1,B2,CN2
720 90 FORMAT(4D25.15)
721 DO 7 I = 1,L
722 DO 7 J = 1,L
723 R1(I,J) = (0.0D0,0.0D0)
724 R2(I,J) = (0.0D0,0.0D0)
725 R3(I,J) = (0.0D0,0.0D0)
726 R4(I,J) = (0.0D0,0.0D0)
727 7 CONTINUE
728 NR = L/2
729 DO 200 I = 1,NR
730 DO 200 J = 1,NR
731 IF(N-I)8,8,200
732 8 IF(N-J)9,9,200
733 9 L1 = I+J+1
734 W1 = (0.0D0,0.0D0)
735 W2 = (0.0D0,0.0D0)
736 W3 = (0.0D0,0.0D0)
737 W4 = (0.0D0,0.0D0)
738 W5 = (0.0D0,0.0D0)
739 W6 = (0.0D0,0.0D0)
740 W7 = (0.0D0,0.0D0)
741 W8 = (0.0D0,0.0D0)
742 DO 100 L = 1,L1
743 K = L-1
744 IF(K-IABS(I-J))100,10,10
745 10 CONTINUE
746 IF(N)11,11,12
747 11 IF(MOD((I+J+K),2).NE.0) GO TO 100
748 FACT1 = 1.0D0
749 GO TO 13
750 12 FACT1 = DFLOAT((-1)**N)
751 13 J1 = 2*I
752 J2 = 2*J
753 J3 = 2*K
754 M1 = 2*N
755 C Seite 32
756
757
758 F3J1 = TRIXJ(J1,J2,J3,M1,FCT,NRISK,NTRIX)
759 FACT2 = 0.5D0*DFLOAT(2*K+1)*F3J1*
760
761
762 1DSQRT(DFLOAT((2*I+1)*(2*J+1))/DFLOAT(I*(I+1)*J*(J+1)))
763 FACT3 = FACT1*FACT2
764 IF(MOD((J-I+K),2).NE.0) GO TO 27
765 IF(J-I+K)25,23,24
766 23 FACT4 = 1.0D0
767 GO TO 26
768 24 FACT4 = DFLOAT((-1)**((J-I+K)/2))
769 GO TO 26
770 25 FACT4 = DFLOAT((-1)**((I-J-K)/2))
771 26 J1 = 2*I
772 J2 = 2*J
773 J3 = 2*K
774 M1 = 0
775 F3J2 = TRIXJ(J1,J2,J3,M1,FCT,NRISK,NTRIX)
776 FACT5 = DFLOAT(I*(I+1)+J*(J+1)-K*(K+1))*F3J2
777 GO TO 28
778 27 FACT6 = 0.0D0
779 GO TO 29
780 28 FACT6 = FACT4*FACT5
781 29 IF(K-IABS(I-J)-1)44,30,30
782 30 IF(MOD((J-I+K+1),2).NE.0) GO TO 44
783 IF(J-I+K+1)33,31,32
784 31 FACT7 = 1.0D0
785 GO TO 41
786 32 FACT7 = DFLOAT((-1)**((J-I+K+1)/2))
787 GO TO 41
788 33 FACT7 = DFLOAT((-1)**((I-J-K-1)/2))
789 41 J1 = 2*I
790 J2 = 2*J
791 J3 = 2*(K-1)
792 M1 = 0
793 F3J3 = TRIXJ(J1,J2,J3,M1,FCT,NRISK,NTRIX)
794 IF(IABS(I-J))42,42,43
795 42 FACT8 = DFLOAT(K)*DSQRT(DFLOAT((I+J-1)**2-K**2))*F3J3
796 GO TO 45
797 43 FACT8 = DSQRT(DFLOAT((K**2-IABS(I-J)**2)*((I+J+1)**2-K**2)))*F3J3
798 GO TO 45
799 44 FACT9 = 0.0D0
800 GO TO 50
801 45 FACT9 = FACT7*FACT8
802 50 CONTINUE
803 IF(K)51,51,52
804 51 FACT10 = 1.0D0
805 GO TO 53
806 52 FACT10 = DFLOAT((-1)**K)
807 53 ACR1 = FACT3*FACT6
808 ACR2 = FACT10*ACR1
809 ACR3 = FACT3*FACT9
810 ACR4 = FACT10*ACR3
811 K1 = K+1
812 W1 = W1+DCMPLX(ACR1*X1(K1),0.0D0)
813 W2 = W2+DCMPLX(ACR1*X1(K1),ACR1*Y1(K1))
814 W3 = W3+DCMPLX(ACR2*X1(K1),ACR2*Y1(K1))
815 W4 = W4+DCMPLX(ACR1*X2(K1),0.0D0)
816 IF(N)100,100,99
817 C Seite 33
818
819
820 99 CONTINUE
821 W5 = W5+DCMPLX(ACR3*X1(K1),0.0D0)
822 W6 = W6+DCMPLX(ACR3*X1(K1),ACR3*Y1(K1))
823 W7 = W7+DCMPLX(ACR4*X1(K1),ACR4*Y1(K1))
824
825
826 W8 = W8+DCMPLX(ACR3*X2(K1),0.0D0)
827 100 CONTINUE
828 R1(2*I-1,2*J-1) = W1
829 R2(2*I-1,2*J-1) = W2
830 R3(2*I-1,2*J-1) = W3
831 R4(2*I-1,2*J-1) = W4
832 R1(2*I,2*J) = W1
833 R2(2*I,2*J) = W2
834 R3(2*I,2*J) = W3
835 R4(2*I,2*J) = W4
836 GO TO (101,103),NP1
837 101 CONTINUE
838 IF(N)200,200,102
839 102 CONTINUE
840 R1(2*I-1,2*J) = W5
841 R2(2*I-1,2*J) = W6
842 R3(2*I-1,2*J) = W7
843 R4(2*I-1,2*J) = W8
844 R1(2*I,2*J-1) = -W5
845 R2(2*I,2*J-1) = -W6
846 R3(2*I,2*J-1) = -W7
847 R4(2*I,2*J-1) = -W8
848 GO TO 200
849 103 CONTINUE
850 IF(N)200,200,104
851 104 CONTINUE
852 R1(2*I-1,2*J) = -W5
853 R2(2*I-1,2*J) = -W6
854 R3(2*I-1,2*J) = -W7
855 R4(2*I-1,2*J) = -W8
856 R1(2*I,2*J-1) = W5
857 R2(2*I,2*J-1) = W6
858 R3(2*I,2*J-1) = W7
859 R4(2*I,2*J-1) = W5
860 200 CONTINUE
861 RETURN
862 END
863 C Seite 34
864
865
866 SUBROUTINE MCNV(A,N,D,L,M)
867 DIMENSION A(1),L(1),M(1)
868 COMPLEX*16 A,D,BIGA,HOLD
869 D = (1.0D0,0.0D0)
870 NK=-N
871 DO 80 K=1,N
872 NK=NK+N
873 L(K)=K
874 M(K)=K
875 KK=NK+K
876 BIGA=A(KK)
877 DO 20 J=K,N
878 IZ=N*(J-1)
879 DO 20 I=K,N
880 IJ=IZ+1
881 10 IF(CDABS(BIGA)-CDABS(A(IJ))) 15,20,20
882 15 BIGA=A(IJ)
883 L(K)=I
884 M(K)=J
885 20 CONTINUE
886 J=L(K)
887 IF(J-K) 35,35,25
888 25 KI=K-N
889 DO 30 I=1,N
890 KI=KI+N
891 HOLD=-A(KI)
892 JI=KI-K+J
893 A(KI)=A(JI)
894 30 A(JI) =HOLD
895 35 I=M(K)
896 IF(I-K) 45,45,38
897 38 JP=N*(I-1)
898 DO 40 J=1,N
899 JK=NK+J
900 JI=JP+J
901 HOLD=-A(JK)
902 A(JK)=A(JI)
903 40 A(JI) =HOLD
904 45 IF(CDABS(BIGA)) 48,46,48
905 46 D = (0.0D0,0.0D0)
906 RETURN
907 48 DO 55 I=1,N
908 IF(I-K) 50,55,50
909 50 IK=NK+1
910 A(IK)=A(IK)/(-BIGA)
911 55 CONTINUE
912 DO 65 I=1,N
913 IK=NK+1
914 HOLD=A(IK)
915 IJ=I-N
916 DO 65 J=1,N
917 IJ=IJ+N
918 IF(I-K) 60,65,60
919 60 IF(J-K) 62,65,62
920 62 KJ=IJ-I+K
921 A(IJ)=HOLD*A(KJ)+A(IJ)
922 65 CONTINUE
923 KJ=K-N
924 DO 75 J=1,N
925 C Seite 35
926
927
928 KJ=KJ+N
929
930
931 IF(J-K) 70,75,70
932 70 A(KJ)=A(KJ)/BIGA
933 75 CONTINUE
934 D=D*BIGA
935 A(KK) = (1.0D0,0.0D0)/BIGA
936 80 CONTINUE
937 K=N
938 100 K=(K-1)
939 IF(K) 150,150,105
940 105 I=L(K)
941
942
943 IF(I-K) 120,120,108
944 108 JQ=N*(K-1)
945 JR=N*(I-1)
946 DO 110 J=1,N
947 JK=JQ+J
948 HOLD=A(JK)
949 JI=JR+J
950 A(JK)=-A(JI)
951 110 A(JI) =HOLD
952 120 J=M(K)
953 IF(J-K) 100,100,125
954 125 KI=K-N
955 DO 130 I=1,N
956 KI=KI+N
957 HOLD=A(KI)
958 JI=KI-K+J
959 A(KI)=-A(JI)
960 130 A(JI) =HOLD
961 GO TO 100
962 150 RETURN
963 END
964 C Seite 36
965
966
967 SUBROUTINE COND(M,ND,RE,IM)
968 DIMENSION RE(ND,1),IM(ND,1)
969 DOUBLE PRECISION RE,IM,SCALE
970 MD = 1
971 IF(M.GT.1) MD = 2*M-1
972 MD1 = MD+1
973 NBGR = ND
974 NROW = NBGR
975 DO 60 KR = MD1,NBGR
976 SCALE = 1.0D0/IM(NROW,NROW)
977 DO 8 LC = MD,NBGR
978 RE(NROW,LC) = SCALE*RE(NROW,LC)
979 IM(NROW,LC) = SCALE*IM(NROW,LC)
980 8 CONTINUE
981 MROW = NROW-1
982 DO 20 MR = MD,MROW
983 SCALE = IM(MR,NROW)
984 DO 16 MC = MD,NBGR
985 RE(MR,MC) = RE(MR,MC)-SCALE*RE(NROW,MC)
986 IM(MR,MC) = IM(MR,MC)-SCALE*IM(NROW,MC)
987 16 CONTINUE
988 20 CONTINUE
989 NROW = NROW-1
990 60 CONTINUE
991 NROW = NBGR-1
992 DO 80 I = 1,NROW
993 IB = I+1
994 DO 72 J = IB,NBGR
995 IM(I,J) = 0.0D0
996 72 CONTINUE
997 80 CONTINUE
998 RETURN
999 END
1000 C Seite 37
1001
1002
1003 SUBROUTINE ORTHO(M,ND,RE,IM,X,Y)
1004 DIMENSION RE(ND,1),IM(ND,1),X(1),Y(1)
1005 DOUBLE PRECISION RE,IM,X,Y,SUM1,SUM2
1006 MD = 1
1007 IF(M.GT.1) MD = 2*M-1
1008 NBGR = ND
1009 SUM1 = 0.0D0
1010 DO 20 K = MD,NBGR
1011 SUM1 = SUM1+RE(NBGR,K)**2+IM(NBGR,K)**2
1012 20 CONTINUE
1013 SUM1 = DSQRT(SUM1)
1014 DO 28 K = MD,NBGR
1015 RE(NBGR,K) = RE(NBGR,K)/SUM1
1016 IM(NBGR,K) = IM(NBGR,K)/SUM1
1017 28 CONTINUE
1018 NMI = NBGR-1
1019 NROW = NBGR
1020 DO 100 I = MD,NMI
1021 NROW = NROW-1
1022 MROW = NROW
1023 DO 36 K = MD,NBGR
1024 X(K) = RE(NROW,K)
1025 Y(K) = IM(NROW,K)
1026 36 CONTINUE
1027 DO 80 J = NROW,NMI
1028 SUM1 = 0.0D0
1029 SUM2 = 0.0D0
1030 MROW = MROW+1
1031 DO 40 K = MD,NBGR
1032 SUM1 = SUM1+RE(MROW,K)*RE(NROW,K)+IM(MROW,K)*IM(NROW,K)
1033 SUM2 = SUM2+RE(MROW,K)*IM(NROW,K)-IM(MROW,K)*RE(NROW,K)
1034 40 CONTINUE
1035 DO 48 K = MD,NBGR
1036 X(K) = X(K)-SUM1*RE(MROW,K)+SUM2*IM(MROW,K)
1037 Y(K) = Y(K)-SUM1*IM(MROW,K)-SUM2*RE(MROW,K)
1038 48 CONTINUE
1039 80 CONTINUE
1040 SUM1 = 0.0D0
1041 DO 84 K = MD,NBGR
1042 SUM1 = SUM1+X(K)**2+Y(K)**2
1043 84 CONTINUE
1044 SUM1 = DSQRT(SUM1)
1045 DO 88 K = MD,NBGR
1046 RE(NROW,K) = X(K)/SUM1
1047 IM(NROW,K) = Y(K)/SUM1
1048 88 CONTINUE
1049 100 CONTINUE
1050 RETURN
1051 END
1052 C Seite 38
1053
1054
1055 SUBROUTINE PERFT(NP,M,NR,ND,AR,AI,BR,BI,T,RE,IM,X,Y)
1056 DIMENSION AR(NR,1),AI(NR,1),BR(NR,1),BI(NR,1),T(ND,1),RE(ND,1),
1057 1IM(ND,1),X(1),Y(1)
1058 DOUBLE PRECISION AR,AI,BR,BI,RE,IM,X,Y,FAC
1059 COMPLEX*16 T
1060 MD = 1
1061 IF(M.GT.1) MD = 2*M-1
1062 NR = ND/2
1063 FAC = 1.0D0
1064 IF(NP.GT.0) FAC = -1.0D0
1065 DO 10 I = 1,NR
1066 DO 10 J = 1,NR
1067 RE(2*I-1,2*J-1) = AR(I,J)
1068 RE(2*I-1,2*J) = FAC*BR(I,J)
1069 RE(2*I,2*J-1) = FAC*BR(I,J)
1070 RE(2*I,2*J) = -AR(I,J)
1071 IM(2*I-1,2*J-1) = AI(I,J)
1072 IM(2*I-1,2*J) = FAC*BI(I,J)
1073 IM(2*I,2*J-1) = FAC*BI(I,J)
1074 IM(2*I,2*J) = -AI(I,J)
1075 IF(1.EQ.J) IM(2*I,2*I) = 1.0D0-AI(I,I)
1076 10 CONTINUE
1077 CALL COND(M,ND,RE,IM)
1078 CALL ORTHO(M,ND,RE,IM,X,Y)
1079 DO 11 I = 1,ND
1080 DO 11 J = 1,ND
1081 T(I,J) = (0.0D0,0.0D0)
1082 11 CONTINUE
1083 DO 12 I = MD,ND
1084 DO 12 J = MD,ND
1085 DO 12 K = MD,ND
1086 T(I,J) = T(I,J)-DCMPLX(RE(K,J)*RE(K,J),-IM(K,I)*RE(K,J))
1087 12 CONTINUE
1088 RETURN
1089 END
1090 C Seite 39
1091
1092
1093 SUBROUTINE LINE(THETA,NIP,C,ALPHA,R,DR)
1094 DOUBLE PRECISION THETA,C,ALPHA,R,DR
1095 R = C*DSIN(ALPHA)/DSIN(THETA+DFLOAT(NIP)*ALPHA)
1096 DR = -R/DTAN(THETA+DFLOAT(NIP)*ALPHA)
1097 RETURN
1098 END
1099
1100
1101 SUBROUTINE TRCIR(STH,CTH,C,A,R,DR)
1102 DOUBLE PRECISION STH,CTH,C,A,R,DR,X,ST,CT
1103 ST = C*STH
1104 CT = C*CTH
1105 X = DSQRT(A**2-ST**2)
1106 R = CT+X
1107 DR = -ST-ST*CT/X
1108 RETURN
1109 END
1110
1111
1112 SUBROUTINE ELLIPS(STH,CTH,AZ,BRA,R,DR)
1113 DOUBLE PRECISION STH,CTH,AZ,BRA,R,DR
1114 R = AZ/DSQRT (CTH**2+(AZ*STH/BRA)**2)
1115 DR = R**3*STH*CTH*(1.0D0-(AZ/BRA)**2)/AZ**2
1116 RETURN
1117 END
1118
1119
1120 SUBROUTINE TRELLI(STH,CTH,AZ,BRA,C,R,DR)
1121 DOUBLE PRECISION STH,CTH,AZ,BRA,C,R,DR,NE,RO
1122 NE = (AZ*STH)**2+(BRA*CTH)**2
1123 RO = NE-(C*STH)**2
1124 RO = DSQRT(RO)
1125 R = (BRA**2*C*CTH+AZ*BRA*RO)/NE
1126 DR = -2.0D0*STH*CTH*(AZ**2-BRA**2)*R/NE+(-BRA**2*C*STH+STH*CTH*AZ*
1127 1BRA*(AZ**2-BRA**2-C**2)/RO)/NE
1128 RETURN
1129 END